##############################
#
# Replication file for the individual analysis in:
#
# The Morning After:
# Cabinet Instability and the Purging of Ministers after Failed Coup Attempts in Autocracies
#
# Published in the Journal of Politics
#
##################

# Set-up

windowsFonts(Times=windowsFont("Times New Roman"))

pacman::p_load(tidyverse,texreg,haven,lfe)

############
# Figure 1 #
############

df_figure1 <- read_dta("morning_after_dataset_countrylevel.dta") %>% 
              filter(democracy1==0 & independent==1 & !is.na(coupattempt_whogov) & !is.na(replacement_rateadj_minister)) %>%
              mutate(coupattempt_whogov = as.factor(coupattempt_whogov),
                     coupattempt_whogov = recode(coupattempt_whogov,"1" = "Failed coup attempt","0" = "No failed coup attempt")) %>% 
              dplyr::select(coupattempt_whogov,replacement_rateadj_minister)
# Plot

figure_1 <- ggplot(df_figure1, aes(x=coupattempt_whogov, y=replacement_rateadj_minister)) + 
  geom_violin(fill='#cccccc',alpha=0.8) + 
  theme_bw() +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),axis.title.y=element_text(size=16,angle=0,vjust=1),
        axis.text=element_text(size=16),
        text=element_text(family="Times")) +
        ylab("Replacement\nrates") +
        xlab("") +
  geom_boxplot(width=0.1)

# Data points

mean_violin <- df_figure1 %>% group_by(coupattempt_whogov) %>% summarize(median = median(replacement_rateadj_minister,na.rm=TRUE),n = n())

###
# Print ---
###

ggsave(
  'output/figure1.jpg',
  figure_1,
  width = 8,
  height = 6,
  dpi = 1200
)

#############################
# Individual level analysis #
#############################

###
## Load data ---
###

df_aut <- read_rds("morning_after_dataset_individuallevel.rds") %>% filter(democracy_lag == 0 & leader == 0 & !is.na(purged_y1) & !is.na(gdp_cap_pwt_ln) & !is.na(pop_pwt_ln))
  
##
# Get Number of coup attempts ---
##

df_aut_n <- df_aut %>% filter(!is.na(coupattempt_presentyear))

df_aut_names <- length(unique(df_aut_n$name))

n_countries <-  length(unique(df_aut_n$country_isocode))

table(df_aut_n$coupattempt_presentyear)

####
## Set up template for graphs ---
####

pd <- position_dodge(0.75)

plot_canvas <- function(df=df,x=name,y=estimate,se=se){
  ggplot(data = df, aes(y = estimate, x = name)) + 
  geom_point(position = pd) + 
  geom_errorbar(aes(ymin = estimate-1.96*se, ymax = estimate+1.96*se),
                width = 0,
                position = pd,
                size = 1) +
  coord_flip() +
  ylim(-0.1,0.5) +
  labs(x = "", 
       y = "Marginal Effects",
       color = "black") +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  theme_bw() +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
        legend.position = "",plot.title = element_text(size=14),
        axis.text.y = element_text(colour = "black",size=14),
        axis.text.x = element_text(colour = "black",size=14),
        text=element_text(family="Times",size=14,color="black"))
}

###
# Figure 3, 4, and 5 ---
###

# Responsibility ---

responsibility_felm <- felm(purged_y1 ~ minister_type/coupattempt_presentyear +
                              gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                              I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                       "Minister of Natural Resources","Other type of minister"),
                              estimate = c(responsibility_felm$coefficients[10],
                                           responsibility_felm$coefficients[11],
                                           responsibility_felm$coefficients[12],
                                           responsibility_felm$coefficients[13],
                                           responsibility_felm$coefficients[14]),
                              se = c(responsibility_felm$coefficients[10,2],
                                     responsibility_felm$coefficients[11,2],
                                     responsibility_felm$coefficients[12,2],
                                     responsibility_felm$coefficients[13,2],
                                     responsibility_felm$coefficients[14,2]
                              )) %>% mutate(name = fct_relevel(name,
                                                               "Minister of Natural Resources",
                                                               "Other type of minister",
                                                               "Minister of Finance",
                                                               "Minister of Foreign Affairs",
                                                               "Minister of Defense"
                              ))


plot_resp <- plot_canvas(coef_responsibility) + ggtitle("H3a: Ministerial responsibility")

# Affiliation ---

df_party <- df_aut %>% filter(party_group != "unknown") %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm <- felm(purged_y1 ~ party_group/coupattempt_presentyear +
                                 gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                                 I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                              estimate = c(affiliation_felm$coefficients[8],
                                           affiliation_felm$coefficients[9],
                                           affiliation_felm$coefficients[10]),
                              se = c(affiliation_felm$coefficients[8,2],
                                     affiliation_felm$coefficients[9,2],
                                     affiliation_felm$coefficients[10,2]
                              )
                              ) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"))

plot_affiliation <- plot_canvas(coef_affiliation) + ggtitle("H2b: Affiliation")

# Importance ---

df_importance <- df_aut %>% filter(importance_factor != "0" & !is.na(importance_factor))

importance_felm <- felm(purged_y1 ~ importance_factor/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance <- tibble(name = c("Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                    "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                           estimate = c(importance_felm$coefficients[10],
                                        importance_felm$coefficients[11],
                                        importance_felm$coefficients[12],
                                        importance_felm$coefficients[13],
                                        importance_felm$coefficients[14]),
                           se = c(importance_felm$coefficients[10,2],
                                  importance_felm$coefficients[11,2],
                                  importance_felm$coefficients[12,2],
                                  importance_felm$coefficients[13,2],
                                  importance_felm$coefficients[14,2]
                           )
                           ) %>% mutate(name = fct_relevel(name,
                                                           "Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                                           "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"))


plot_importance <- plot_canvas(coef_importance) + ggtitle("H3b: Importance in cabinet") + ylim(-0.1,0.6)

# Experience ---

df_experience <- df_aut %>% filter(importance_factor != "0.5")

quantiles_experience <- df_aut %>% summarise(x = quantile(experience, c(0.25, 0.5, 0.75)), experience = c(0.25, 0.5, 0.75))

experience_group <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_presentyear +
                        gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_group <- tibble(name = c("Less than 2 years","2-3 years",
                                    "4-6 years","Over 6 years"),
                           estimate = c(experience_group$coefficients[6],
                                        experience_group$coefficients[7],
                                        experience_group$coefficients[8],
                                        experience_group$coefficients[9]),
                           se = c(experience_group$coefficients[6,2],
                                  experience_group$coefficients[7,2],
                                  experience_group$coefficients[8,2],
                                  experience_group$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "Over 6 years","4-6 years",
                                                            "2-3 years",
                                                            "Less than 2 years"
                                                            
                           ))


plot_expgroup <- plot_canvas(coef_experience_group) + ggtitle("H2a: Experience")

# Combination (typology) ---

medianbyyear <- df_experience %>% group_by(country_isocode,year) %>% summarize(experience_median = median(experience,na.rm=TRUE))

df_combination <- df_experience %>% left_join(.,medianbyyear,by=c("country_isocode","year")) %>% filter(party_group != "unknown") %>%
                                                                                                 mutate(experience_bin = ifelse(experience >= experience_median, 1,0),
                                                                                                 important_bin = if_else(importance_numeric %in% c(3,4),1,0),
                                                                                                 loyal_bin = if_else(experience_bin == 1 & party_group == "leader_party",1,0),
                                                                                                 typology = case_when(important_bin == 0 & loyal_bin == 0 ~ 1,
                                                                                                                      important_bin == 1 & loyal_bin == 0 ~ 2,
                                                                                                                      important_bin == 0 & loyal_bin == 1 ~ 3,
                                                                                                                      important_bin == 1 & loyal_bin == 1 ~ 4
                                                                                                 )
)

combination_felm <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                    "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                              estimate = c(combination_felm$coefficients[6],
                                           combination_felm$coefficients[7],
                                           combination_felm$coefficients[8],
                                           combination_felm$coefficients[9]),
                              se = c(combination_felm$coefficients[6,2],
                                     combination_felm$coefficients[7,2],
                                     combination_felm$coefficients[8,2],
                                     combination_felm$coefficients[9,2]
                              )) %>% mutate(name = fct_relevel(name,
                                                               "High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                               "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"
                                                               
                              ))


plot_combination <- plot_canvas(coef_combination) + ggtitle("H4: Combinations of traits")

# Get N's

n_exp <-  df_experience %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
                             mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_presentyear) & !is.na(minister_type)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

# Print ---

ggsave(
  'output/figure3.jpg',
  gridExtra::grid.arrange(plot_expgroup, plot_affiliation,ncol=2),
  width = 13,
  height = 6,
  dpi = 120
)

ggsave(
  'output/figure4.jpg',
  gridExtra::grid.arrange(plot_resp, plot_importance,ncol=2),
  width = 13,
  height = 6,
  dpi = 120
)

ggsave(
  'output/figure5.jpg',
  plot_combination,
  width = 10,
  height = 6,
  dpi = 120
)

################
## Appendix K ##
################

# Responsibility ---

df_aut <- df_aut %>% mutate(minister_type = fct_relevel(minister_type,
                                                      "Minister of Defense",
                                                      "Minister of Foreign Affairs",
                                                      "Minister of Finance",
                                                      "Other type of minister",                                                      
                                                      "Minister of Natural Resources"
                                                      ))


responsibility_felm <- felm(purged_y1 ~ minister_type*coupattempt_presentyear +
                              gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                              I(experience^3) | country_isocode+year | 0 | country_isocode, data=df_aut)


# Party affiliation ---

df_party <- df_party %>% mutate(party_group = fct_relevel(party_group,
                                                          "From another party",
                                                          "No party affiliation",
                                                            "From the leader's party"
))

affiliation_felm <- felm(purged_y1 ~ party_group*coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode, data=df_party)

# Importance ---

df_importance <- df_importance %>% mutate(importance_factor = fct_relevel(importance_factor,
                                                             "Prime minister/President\n(not leader)",
                                                             "Vice president,\ndeputy prime minister,\ntop minister",
                                                             "Medium-ranking\nminister",
                                                             "Low-ranking\nminister",
                                                             "Junior minister\nor other low-\nranking post"
))

importance_felm <- felm(purged_y1 ~ importance_factor*coupattempt_presentyear +
                          gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                          I(experience^3) | country_isocode+year | 0 | country_isocode, data=df_importance)

# Experience ---

experience_felm <- felm(purged_y1 ~ experience_group*coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience)


# Combination ---

combination_felm <- felm(purged_y1 ~ relevel(as.factor(typology),2)*coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination)


# Tables ---

# Table for Figure 3

texreg(list(experience_felm,affiliation_felm),stars = c(0.05),
       file="output/appendix_k1.tex",
       custom.model.names=c("H2a: Experience","H2b: Affiliation"),
       custom.coef.names=c("2-3 years of experience (Ref: <2 y)","4-6 years of experience (Ref: <2 y)","Over 6 years of experience (Ref: <2 y)",
                           "Failed coup attempt", "Log of GDP per capita","Log of Population",
                           "FCA*2-3 years of experience (Ref: <2 y)","FCA*4-6 years of experience (Ref: <2 y)","FCA*Over 6 years of experience (Ref: <2 y)",
                           "No party affiliation (Ref: Other party)","From the leader's party (Ref: Other party)","Experience","Experience^2","Experience^3",
                           "FCA*No party affiliation (Ref: Other party)",
                           "FCA*From the leader's party (Ref: Other party)"),
       reorder.coef=c(1,2,3,10,11,4,7,8,9,15,16,5,6,12,13,14),
       label = "indtable",
       single.row = TRUE,
       caption ="Table for Figure 3",
       caption.above=TRUE,
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE,
       custom.note = ("\\parbox{\\linewidth}{\\vspace{2pt}%stars. Dependent variable: Purged next year. 
                      All models include country and year fixed effects. Country-clustered standard errors in parentheses. 
                      FCA = Failed coup attempt}"))

# Table for Figure 4

texreg(list(responsibility_felm,importance_felm),stars = c(0.05),
       file="output/appendix_k2.tex",
       custom.model.names=c("H3a: Responsibility","H3b: Importance"),
       custom.coef.names=c("Minister of Foreign Affairs (Ref: D)","Minister of Finance (Ref: D)", "Other type of minister (Ref: D)",
                            "Minister of Natural Resources (Ref: D)","Failed coup attempt","Log of GDP per capita","Log of Population",
                            "Experience","Experience^2","Experience^3","FCA: Minister of Foreign Affairs (Ref: D)",
                            "FCA*Minister of Finance (Ref: D)","FCA*Other type of minister (Ref: D)","FCA*Minister of Natural Resources (Ref: D)",
                            "VP, DP, top minister (Ref: PMP)","Medium-ranking minister (Ref: PMP)","Low-ranking minister (Ref: PMP)",
                           "Junior minister (Ref: PMP)","FCA*VP, DP, top minister (Ref: PMP)",
                           "FCA*Medium-ranking minister (Ref: PMP)","FCA*Low-ranking minister (Ref: PMP)",
                           "FCA*Junior minister (Ref: PMP)"),
       reorder.coef=c(1,2,3,4,
                      15,16,17,18,
                      5,
                      11,12,13,14,
                      19,20,21,22,
                      6,7,8,9,10),
       label = "indtable2",
       single.row = TRUE,
       caption ="Table for Figure 4",
       caption.above=TRUE,
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE,
       custom.note = ("\\parbox{\\linewidth}{\\vspace{2pt}%stars. Dependent variable: Purged next year. 
                      All models include country and year fixed effects. Country-clustered standard errors in parentheses.
                      FCA = Failed coup attempt, D = Defense, PMP = Prime minister or President.}"))

# Table for Figure 5

texreg(list(combination_felm),stars = c(0.05),
       file="output/appendix_k3.tex",
       custom.model.names=c("Combination of traits"),
       custom.coef.names=c("Low responsibility and weak signs of loyalty (Ref: HW)",
                           "Low responsibility and strong signs of loyalty (Ref: HW)",
                          "High responsibility and strong signs of loyalty (Ref: HW)",
                          "Failed coup attempt","Log of GDP per capita","Log of Population",
                          "FCA*High responsibility and weak signs of loyalt (Ref: HW)",
                          "FCA*Low responsibility and strong signs of loyalty (Ref: HW)",
                          "FCA*High responsibility and strong signs of loyalty (Ref: HW)"),
       reorder.coef=c(1,2,3,4,
                      7,8,9,
                      6,5),
       label = "indtable3",
       single.row = TRUE,
       caption ="Table for Figure 5",
       caption.above=TRUE,
       include.ci = FALSE,
       include.adjrs = FALSE,
       include.rmse = FALSE,
       custom.note = ("\\parbox{\\linewidth}{\\vspace{2pt}%stars. Dependent variable: Purged next year. 
                      All models include country and year fixed effects. Country-clustered standard errors in parentheses.
                      FCA = Failed coup attempt, LW = High responsibility and weak signs of loyalty.}"))

################
## Appendix O ##
################

# Affiliation ---

df_party_morethanone <- df_aut %>% filter(party_group != "unknown" & n_party > 0) %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm_morethanone <- felm(purged_y1 ~ party_group/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party_morethanone) %>% summary()

coef_affiliation_morethanone <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(affiliation_felm_morethanone$coefficients[8],
                                        affiliation_felm_morethanone$coefficients[9],
                                        affiliation_felm_morethanone$coefficients[10]),
                           se = c(affiliation_felm_morethanone$coefficients[8,2],
                                  affiliation_felm_morethanone$coefficients[9,2],
                                  affiliation_felm_morethanone$coefficients[10,2]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"))

plot_affiliation_morethanon <- plot_canvas(coef_affiliation_morethanone) + ggtitle("H2b: Affiliation")

# Combination (typology) ---

df_combination_morethanone <- df_combination %>% filter(party_group != "unknown" & n_party > 0)

combination_felm_morethanone <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination_morethanone) %>% summary()

coef_combination_morethanone <- tibble(name = c("Low responsibility and\nweak signs of loyalty","High responsibility and\nweak signs of loyalty",
                                    "Low responsibility and\nstrong signs of loyalty","High responsibility and\nstrong signs of loyalty"),
                           estimate = c(combination_felm_morethanone$coefficients[6],
                                        combination_felm_morethanone$coefficients[7],
                                        combination_felm_morethanone$coefficients[8],
                                        combination_felm_morethanone$coefficients[9]),
                           se = c(combination_felm_morethanone$coefficients[6,2],
                                  combination_felm_morethanone$coefficients[7,2],
                                  combination_felm_morethanone$coefficients[8,2],
                                  combination_felm_morethanone$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "High responsibility and\nstrong signs of loyalty","Low responsibility and\nstrong signs of loyalty",
                                                            "High responsibility and\nweak signs of loyalty","Low responsibility and\nweak signs of loyalty"
                                                            
                           ))


plot_combination_morethanone <- plot_canvas(coef_combination_morethanone) + ggtitle("H4: Combinations of traits")

# Get N's ---

n_exp <-  df_party_morethanone %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination_morethanone %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

###
# Print ---
###

ggsave(
  'output/appendix_o1.jpg',
  gridExtra::grid.arrange(plot_affiliation_morethanon, plot_combination_morethanone,ncol=2),
  width = 10,
  height = 6,
  dpi = 120
)

################
## Appendix L ##
################

# Responsiblity ---

responsibility_felm_br <- felm(purged_y1 ~ minister_type/coupattempt_presentyear_br +
                              gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                              I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility_br <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                       "Minister of Natural Resources","Other type of minister"),
                              estimate = c(responsibility_felm_br$coefficients[10],
                                           responsibility_felm_br$coefficients[11],
                                           responsibility_felm_br$coefficients[12],
                                           responsibility_felm_br$coefficients[13],
                                           responsibility_felm_br$coefficients[14]),
                              se = c(responsibility_felm_br$coefficients[10,2],
                                     responsibility_felm_br$coefficients[11,2],
                                     responsibility_felm_br$coefficients[12,2],
                                     responsibility_felm_br$coefficients[13,2],
                                     responsibility_felm_br$coefficients[14,2]
                              )) %>% mutate(name = fct_relevel(name,
                                                               "Minister of Natural Resources",
                                                               "Other type of minister",
                                                               "Minister of Finance",
                                                               "Minister of Foreign Affairs",
                                                               "Minister of Defense"
                              ))


plot_resp_br <- plot_canvas(coef_responsibility_br) + ggtitle("H3a: Ministerial responsibility")

# Affiliation ---

df_party <- df_aut %>% filter(party_group != "unknown") %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm_br <- felm(purged_y1 ~ party_group/coupattempt_presentyear_br +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation_br <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(affiliation_felm_br$coefficients[8],
                                        affiliation_felm_br$coefficients[9],
                                        affiliation_felm_br$coefficients[10]),
                           se = c(affiliation_felm_br$coefficients[8,2],
                                  affiliation_felm_br$coefficients[9,2],
                                  affiliation_felm_br$coefficients[10,2]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"))

plot_affiliation_br <- plot_canvas(coef_affiliation_br) + ggtitle("H2b: Affiliation")

# Importance ---

df_importance <- df_aut %>% filter(importance_factor != "0")

importance_felm_br <- felm(purged_y1 ~ importance_factor/coupattempt_presentyear_br +
                          gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                          I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance_br <- tibble(name = c("Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                   "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                          estimate = c(importance_felm_br$coefficients[10],
                                       importance_felm_br$coefficients[11],
                                       importance_felm_br$coefficients[12],
                                       importance_felm_br$coefficients[13],
                                       importance_felm_br$coefficients[14]),
                          se = c(importance_felm_br$coefficients[10,2],
                                 importance_felm_br$coefficients[11,2],
                                 importance_felm_br$coefficients[12,2],
                                 importance_felm_br$coefficients[13,2],
                                 importance_felm_br$coefficients[14,2]
                          )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"))


plot_importance_br <- plot_canvas(coef_importance_br) + ggtitle("H3b: Importance in cabinet") + ylim(-0.1,0.6)

# Experience ---

experience_group_br <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_presentyear_br +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_group <- tibble(name = c("Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years"),
                                estimate = c(experience_group_br$coefficients[6],
                                             experience_group_br$coefficients[7],
                                             experience_group_br$coefficients[8],
                                             experience_group_br$coefficients[9]),
                                se = c(experience_group_br$coefficients[6,2],
                                       experience_group_br$coefficients[7,2],
                                       experience_group_br$coefficients[8,2],
                                       experience_group_br$coefficients[9,2]
                                )) %>% mutate(name = fct_relevel(name,
                                                                 "Over 6 years","4-6 years",
                                                                 "2-3 years",
                                                                 "Less than 2 years"
                                                                 
                                ))


plot_expgroup_br <- plot_canvas(coef_experience_group) + ggtitle("H2a: Experience")

# Combination ---

combination_felm_br <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear_br +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination_br <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                    "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                           estimate = c(combination_felm_br$coefficients[6],
                                        combination_felm_br$coefficients[7],
                                        combination_felm_br$coefficients[8],
                                        combination_felm_br$coefficients[9]),
                           se = c(combination_felm_br$coefficients[6,2],
                                  combination_felm_br$coefficients[7,2],
                                  combination_felm_br$coefficients[8,2],
                                  combination_felm_br$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                            "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"
                                                            
                           ))


plot_combination_br <- plot_canvas(coef_combination_br) + ggtitle("H4: Combinations of traits")

# Get N's ---

n_exp <-  df_experience %>% filter(!is.na(coupattempt_presentyear_br)) %>% group_by(coupattempt_presentyear_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_presentyear_br)) %>% group_by(coupattempt_presentyear_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_presentyear_br) & !is.na(minister_type)) %>% group_by(coupattempt_presentyear_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_presentyear_br)) %>% group_by(coupattempt_presentyear_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_presentyear_br)) %>% group_by(coupattempt_presentyear_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

# Print ---

ggsave(
  'output/appendix_l1.jpg',
  gridExtra::grid.arrange(plot_expgroup_br, plot_affiliation_br, plot_resp_br, plot_importance_br),
  width = 10,
  height = 9,
  dpi = 120
)

ggsave(
  'output/appendix_l2.jpg',
  plot_combination_br,
  width = 8,
  height = 9,
  dpi = 120
)

################
## Appendix P ##
################

## Set up template ---

pd <- position_dodge(0.75)

plot_canvas_type <- function(df=df,x=name,y=estimate,se=se, group = group){
    ggplot(data = df, aes(y = estimate, x = name, group = group, colour = group)) + 
    geom_point(position = pd) + 
    geom_errorbar(aes(ymin = estimate-1.96*se, ymax = estimate+1.96*se),
                  width = 0,
                  position = pd,
                  size = 1) +
    coord_flip() +
    ylim(-0.2,0.5) +
    labs(x = "", 
         y = "Marginal Effects",
         color = "black") +
    geom_hline(yintercept = 0,
               linetype = "dashed") +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
          legend.position = "bottom",plot.title = element_text(size=14),
          axis.text.y = element_text(colour = "black",size=14),
          axis.text.x = element_text(colour = "black",size=14),
          text=element_text(family="Times",size=14,color="black")) +
    scale_colour_manual(
      name = "", 
      breaks = c("Civilian coup","Military coup"),
      labels = c("Civilian coup","Military coup"),
      values = c("Grey", "Black"))
}

###
# Marginal effects ---
###

# Responsiblity ---

responsibility_felm_brtype <- felm(purged_y1 ~ minister_type/coupattempt_whogov_type_br +
                                 gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                                 I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility_br_type <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                          "Minister of Natural Resources","Other type of minister",
                                          "Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                          "Minister of Natural Resources","Other type of minister"),
                                 estimate = c(responsibility_felm_brtype$coefficients[10],
                                              responsibility_felm_brtype$coefficients[11],
                                              responsibility_felm_brtype$coefficients[12],
                                              responsibility_felm_brtype$coefficients[13],
                                              responsibility_felm_brtype$coefficients[14],
                                              responsibility_felm_brtype$coefficients[15],
                                              responsibility_felm_brtype$coefficients[16],
                                              responsibility_felm_brtype$coefficients[17],
                                              responsibility_felm_brtype$coefficients[18],
                                              responsibility_felm_brtype$coefficients[19]),
                                 se = c(responsibility_felm_brtype$coefficients[10,2],
                                        responsibility_felm_brtype$coefficients[11,2],
                                        responsibility_felm_brtype$coefficients[12,2],
                                        responsibility_felm_brtype$coefficients[13,2],
                                        responsibility_felm_brtype$coefficients[14,2],
                                        responsibility_felm_brtype$coefficients[15,2],
                                        responsibility_felm_brtype$coefficients[16,2],
                                        responsibility_felm_brtype$coefficients[17,2],
                                        responsibility_felm_brtype$coefficients[18,2],
                                        responsibility_felm_brtype$coefficients[19,2]),
                                 group = c(rep("Civilian coup",5),rep("Military coup",5))) %>% 
                                 mutate(name = fct_relevel(name,
                                                                  "Minister of Natural Resources",
                                                                  "Other type of minister",
                                                                  "Minister of Finance",
                                                                  "Minister of Foreign Affairs",
                                                                  "Minister of Defense"
                                 ),
                                 group = fct_relevel(group,"Military coup","Civilian coup")
                                 )


plot_resp_br_type <- plot_canvas_type(coef_responsibility_br_type) + ggtitle("H3a: Ministerial responsibility")

# Affiliation ---

affiliation_felm_br <- felm(purged_y1 ~ party_group/coupattempt_whogov_type_br +
                              gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                              I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation_br_type <- tibble(name = c("From another party","From the leader's party","No party affiliation",
                                       "From another party","From the leader's party","No party affiliation"),
                              estimate = c(affiliation_felm_br$coefficients[8],
                                           affiliation_felm_br$coefficients[9],
                                           affiliation_felm_br$coefficients[10],
                                           affiliation_felm_br$coefficients[11],
                                           affiliation_felm_br$coefficients[12],
                                           affiliation_felm_br$coefficients[13]),
                              se = c(affiliation_felm_br$coefficients[8,2],
                                     affiliation_felm_br$coefficients[9,2],
                                     affiliation_felm_br$coefficients[10,2],
                                     affiliation_felm_br$coefficients[11,2],
                                     affiliation_felm_br$coefficients[12,2],
                                     affiliation_felm_br$coefficients[13,2]),
                                     group = c(rep("Civilian coup",3),rep("Military coup",3)
                              )) %>% 
              mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"
                                ),
             group = fct_relevel(group,"Military coup","Civilian coup")
             )

plot_affiliation_br_type <- plot_canvas_type(coef_affiliation_br_type) + ggtitle("H2b: Affiliation")

# Importance ---

importance_felm_br <- felm(purged_y1 ~ importance_factor/coupattempt_whogov_type_br +
                             gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                             I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance_br_type <- tibble(name = c("Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                      "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)",
                                      "Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                      "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                             estimate = c(importance_felm_br$coefficients[10],
                                          importance_felm_br$coefficients[11],
                                          importance_felm_br$coefficients[12],
                                          importance_felm_br$coefficients[13],
                                          importance_felm_br$coefficients[14],
                                          importance_felm_br$coefficients[15],
                                          importance_felm_br$coefficients[16],
                                          importance_felm_br$coefficients[17],
                                          importance_felm_br$coefficients[18],
                                          importance_felm_br$coefficients[19]),
                             se = c(importance_felm_br$coefficients[10,2],
                                    importance_felm_br$coefficients[11,2],
                                    importance_felm_br$coefficients[12,2],
                                    importance_felm_br$coefficients[13,2],
                                    importance_felm_br$coefficients[14,2],
                                    importance_felm_br$coefficients[15,2],
                                    importance_felm_br$coefficients[16,2],
                                    importance_felm_br$coefficients[17,2],
                                    importance_felm_br$coefficients[18,2],
                                    importance_felm_br$coefficients[19,2]),
                             group = c(rep("Civilian coup",5),rep("Military coup",5)
                             )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
             group = fct_relevel(group,"Military coup","Civilian coup"))


plot_importance_br_type <- plot_canvas_type(coef_importance_br_type) + ggtitle("H3b: Importance in cabinet") + ylim(-0.2,0.6)

# Experience ---

experience_group_br <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_whogov_type_br +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_group_type <- tibble(name = c("Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years",
                                         "Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years"),
                                estimate = c(experience_group_br$coefficients[6],
                                             experience_group_br$coefficients[7],
                                             experience_group_br$coefficients[8],
                                             experience_group_br$coefficients[9],
                                             experience_group_br$coefficients[10],
                                             experience_group_br$coefficients[11],
                                             experience_group_br$coefficients[12],
                                             experience_group_br$coefficients[13]),
                                se = c(experience_group_br$coefficients[6,2],
                                       experience_group_br$coefficients[7,2],
                                       experience_group_br$coefficients[8,2],
                                       experience_group_br$coefficients[9,2],
                                       experience_group_br$coefficients[10,2],
                                       experience_group_br$coefficients[11,2],
                                       experience_group_br$coefficients[12,2],
                                       experience_group_br$coefficients[13,2]),
                                group = c(rep("Civilian coup",4),rep("Military coup",4))) %>% 
                                mutate(name = fct_relevel(name,
                                                                 "Over 6 years","4-6 years",
                                                                 "2-3 years",
                                                                 "Less than 2 years"
                                                                 
                                ),
                                group = fct_relevel(group,"Military coup","Civilian coup"))


plot_expgroup_br_type <- plot_canvas_type(coef_experience_group_type) + ggtitle("H2a: Experience")

# Combination ---

combination_felm_br <- felm(purged_y1 ~ as.factor(typology)/coupattempt_whogov_type_br +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination_br_type <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                       "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty",
                                       "Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                       "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                              estimate = c(combination_felm_br$coefficients[6],
                                           combination_felm_br$coefficients[7],
                                           combination_felm_br$coefficients[8],
                                           combination_felm_br$coefficients[9],
                                           combination_felm_br$coefficients[10],
                                           combination_felm_br$coefficients[11],
                                           combination_felm_br$coefficients[12],
                                           combination_felm_br$coefficients[13]),
                              se = c(combination_felm_br$coefficients[6,2],
                                     combination_felm_br$coefficients[7,2],
                                     combination_felm_br$coefficients[8,2],
                                     combination_felm_br$coefficients[9,2],
                                     combination_felm_br$coefficients[10,2],
                                     combination_felm_br$coefficients[11,2],
                                     combination_felm_br$coefficients[12,2],
                                     combination_felm_br$coefficients[13,2]),
                              group = c(rep("Civilian coup",4),rep("Military coup",4))) %>% 
                              mutate(name = fct_relevel(name,"High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                             "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"),
                              group = fct_relevel(group,"Military coup","Civilian coup"))

plot_combination_br_type <- plot_canvas_type(coef_combination_br_type) + ggtitle("H4: Combinations of traits")

# Get N's ---

n_exp <-  df_experience %>% filter(!is.na(coupattempt_whogov_type_br) & coupattempt_whogov_type_br !="Royal") %>% group_by(coupattempt_whogov_type_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_whogov_type_br) & coupattempt_whogov_type_br !="Royal") %>% group_by(coupattempt_whogov_type_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_whogov_type_br) & coupattempt_whogov_type_br !="Royal" & !is.na(minister_type)) %>% group_by(coupattempt_whogov_type_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_whogov_type_br) & coupattempt_whogov_type_br !="Royal") %>% group_by(coupattempt_whogov_type_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_whogov_type_br) & coupattempt_whogov_type_br !="Royal") %>% group_by(coupattempt_whogov_type_br) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

###
# Print ---
###

ggsave(
  'output/appendix_p1.jpg',
  gridExtra::grid.arrange(plot_expgroup_br_type, plot_affiliation_br_type, plot_resp_br_type, plot_importance_br_type),
  width = 10,
  height = 9,
  dpi = 120
)

ggsave(
  'output/appendix_p2.jpg',
  plot_combination_br_type,
  width = 8,
  height = 9,
  dpi = 120
)

################
## Appendix R ##
################

###
# Marginal effects ---
###

# Gender

df_aut_gender <- df_aut %>% filter(gender != "Unknown")

hired_gender_felm <- felm(hired ~ gender/coup_lastyear +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_aut_gender) %>% summary()

hired_gender <- tibble(name = c("Female","Male"),
                           estimate = c(hired_gender_felm$coefficients[4],
                                        hired_gender_felm$coefficients[5]),
                           se = c(hired_gender_felm$coefficients[4],
                                  hired_gender_felm$coefficients[5]
                           )
) %>% mutate(name = fct_relevel(name,
                                "Female",
                                "Male"))

hired_plot_gender <- plot_canvas(hired_gender) + ggtitle("Gender") + ylim(-0.2,0.5)

# Party affiliation

hired_affiliation_felm <- felm(hired ~ party_group/coup_lastyear +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

hired_coef_affiliation <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(hired_affiliation_felm$coefficients[5],
                                        hired_affiliation_felm$coefficients[6],
                                        hired_affiliation_felm$coefficients[7]),
                           se = c(hired_affiliation_felm$coefficients[5],
                                  hired_affiliation_felm$coefficients[6],
                                  hired_affiliation_felm$coefficients[7]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"))

hired_plot_affiliation <- plot_canvas(hired_coef_affiliation) + ggtitle("Affiliation") + ylim(-0.3,0.8)

# Position

responsibility_felm_hired <- felm(hired ~ minister_type/coup_lastyear +
                              gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility_hired <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                          "Minister of Natural Resources","Other type of minister"),
                                 estimate = c(responsibility_felm_hired$coefficients[7],
                                              responsibility_felm_hired$coefficients[8],
                                              responsibility_felm_hired$coefficients[9],
                                              responsibility_felm_hired$coefficients[10],
                                              responsibility_felm_hired$coefficients[11]),
                                 se = c(responsibility_felm_hired$coefficients[7],
                                        responsibility_felm_hired$coefficients[8],
                                        responsibility_felm_hired$coefficients[9],
                                        responsibility_felm_hired$coefficients[10],
                                        responsibility_felm_hired$coefficients[11]
                                 )) %>% mutate(name = fct_relevel(name,
                                                                  "Minister of Natural Resources",
                                                                  "Other type of minister",
                                                                  "Minister of Finance",
                                                                  "Minister of Foreign Affairs",
                                                                  "Minister of Defense"
                                 ))


hired_plot_resp <- plot_canvas(coef_responsibility_hired) + ggtitle("Ministerial responsibility") + ylim(-0.2,0.8)

# Importance

df_importance <- df_aut %>% filter(importance_factor != "0")

importance_felm_hired <- felm(hired ~ importance_factor/coup_lastyear +
                          gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance_hired <- tibble(name = c("Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                   "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                          estimate = c(importance_felm_hired$coefficients[7],
                                       importance_felm_hired$coefficients[8],
                                       importance_felm_hired$coefficients[9],
                                       importance_felm_hired$coefficients[10],
                                       importance_felm_hired$coefficients[11]),
                          se = c(importance_felm_hired$coefficients[7],
                                 importance_felm_hired$coefficients[8],
                                 importance_felm_hired$coefficients[9],
                                 importance_felm_hired$coefficients[10],
                                 importance_felm_hired$coefficients[11]
                          )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"))


hired_plot_importance <- plot_canvas(coef_importance_hired) + ggtitle("Importance in cabinet") + ylim(-0.2,0.6)

# Get N's ---

n_gender <-  df_aut_gender %>% filter(!is.na(coup_lastyear) & !is.na(gender)) %>% group_by(coup_lastyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_party <-  df_party %>% filter(!is.na(coup_lastyear) & !is.na(party_group)) %>% group_by(coup_lastyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_responsibility <-  df_aut %>% filter(!is.na(coup_lastyear) & !is.na(minister_type)) %>% group_by(coup_lastyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coup_lastyear)) %>% group_by(coup_lastyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

###
# Print ---
###

ggsave(
  'output/appendix_r1.jpg',
  gridExtra::grid.arrange(hired_plot_gender, hired_plot_affiliation,hired_plot_resp,hired_plot_importance,ncol=2),
  width = 10,
  height = 9,
  dpi = 120
)


################
## Appendix N ##
################

df_aut <- read_rds("morning_after_dataset_individuallevel.rds") %>% 
          filter(democracy_lag == 0 & leader == 0 & !is.na(purged_y1) & !is.na(gdp_cap_pwt_ln) & !is.na(pop_pwt_ln) & leaderexperience_continuous > 1)

###
# Marginal effects ---
###

# Responsibility ---

responsibility_felm <- felm(purged_y1 ~ minister_type/coupattempt_presentyear +
                              gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                              I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                       "Minister of Natural Resources","Other type of minister"),
                              estimate = c(responsibility_felm$coefficients[10],
                                           responsibility_felm$coefficients[11],
                                           responsibility_felm$coefficients[12],
                                           responsibility_felm$coefficients[13],
                                           responsibility_felm$coefficients[14]),
                              se = c(responsibility_felm$coefficients[10,2],
                                     responsibility_felm$coefficients[11,2],
                                     responsibility_felm$coefficients[12,2],
                                     responsibility_felm$coefficients[13,2],
                                     responsibility_felm$coefficients[14,2]
                              )) %>% mutate(name = fct_relevel(name,
                                                               "Minister of Natural Resources",
                                                               "Other type of minister",
                                                               "Minister of Finance",
                                                               "Minister of Foreign Affairs",
                                                               "Minister of Defense"
                              ))


plot_resp <- plot_canvas(coef_responsibility) + ggtitle("H3a: Ministerial responsibility") + ylim(-0.2,0.5)

# Affiliation ---

df_party <- df_aut %>% filter(party_group != "unknown") %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm <- felm(purged_y1 ~ party_group/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(affiliation_felm$coefficients[8],
                                        affiliation_felm$coefficients[9],
                                        affiliation_felm$coefficients[10]),
                           se = c(affiliation_felm$coefficients[8,2],
                                  affiliation_felm$coefficients[9,2],
                                  affiliation_felm$coefficients[10,2]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"))

plot_affiliation <- plot_canvas(coef_affiliation) + ggtitle("H2b: Affiliation") + ylim(-0.2,0.5)

# Importance ---

df_importance <- df_aut %>% filter(importance_factor != "0")

importance_felm <- felm(purged_y1 ~ importance_factor/coupattempt_presentyear +
                          gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                          I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance <- tibble(name = c("Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                   "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                          estimate = c(importance_felm$coefficients[10],
                                       importance_felm$coefficients[11],
                                       importance_felm$coefficients[12],
                                       importance_felm$coefficients[13],
                                       importance_felm$coefficients[14]),
                          se = c(importance_felm$coefficients[10,2],
                                 importance_felm$coefficients[11,2],
                                 importance_felm$coefficients[12,2],
                                 importance_felm$coefficients[13,2],
                                 importance_felm$coefficients[14,2]
                          )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"))


plot_importance <- plot_canvas(coef_importance) + ggtitle("H3b: Importance in cabinet") + ylim(-0.2,0.7)

# Experience ---

df_experience <- df_aut %>% filter(importance_factor != "0.5")

quantiles_experience <- df_aut %>% summarise(x = quantile(experience, c(0.25, 0.5, 0.75)), experience = c(0.25, 0.5, 0.75))

experience_group <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_group <- tibble(name = c("Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years"),
                                estimate = c(experience_group$coefficients[6],
                                             experience_group$coefficients[7],
                                             experience_group$coefficients[8],
                                             experience_group$coefficients[9]),
                                se = c(experience_group$coefficients[6,2],
                                       experience_group$coefficients[7,2],
                                       experience_group$coefficients[8,2],
                                       experience_group$coefficients[9,2]
                                )) %>% mutate(name = fct_relevel(name,
                                                                 "Over 6 years","4-6 years",
                                                                 "2-3 years",
                                                                 "Less than 2 years"
                                                                 
                                ))


plot_expgroup <- plot_canvas(coef_experience_group) + ggtitle("H2a: Experience")

# Combination ---

medianbyyear <- df_experience %>% group_by(country_isocode,year) %>% summarize(experience_median = median(experience,na.rm=TRUE))

df_combination <- df_experience %>% left_join(.,medianbyyear,by=c("country_isocode","year")) %>% filter(party_group != "unknown") %>%
                                                                                                        mutate(experience_bin = ifelse(experience >= experience_median, 1,0),
                                                                                                        important_bin = if_else(importance_numeric %in% c(3,4),1,0),
                                                                                                        loyal_bin = if_else(experience_bin == 1 & party_group == "leader_party",1,0),
                                                                                                        typology = case_when(important_bin == 0 & loyal_bin == 0 ~ 1,
                                                                                                                             important_bin == 1 & loyal_bin == 0 ~ 2,
                                                                                                                             important_bin == 0 & loyal_bin == 1 ~ 3,
                                                                                                                             important_bin == 1 & loyal_bin == 1 ~ 4
                                                                                                        )
)

combination_felm <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                    "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                           estimate = c(combination_felm$coefficients[6],
                                        combination_felm$coefficients[7],
                                        combination_felm$coefficients[8],
                                        combination_felm$coefficients[9]),
                           se = c(combination_felm$coefficients[6,2],
                                  combination_felm$coefficients[7,2],
                                  combination_felm$coefficients[8,2],
                                  combination_felm$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                            "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"
                                                            
                           ))


plot_combination <- plot_canvas(coef_combination) + ggtitle("H4: Combinations of traits") + ylim(-0.2,0.5)

# Get N's ---

n_exp <-  df_experience %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_presentyear) & !is.na(minister_type)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

###
# Print ---
###

ggsave(
  'output/appendix_n1.jpg',
  gridExtra::grid.arrange(plot_expgroup, plot_affiliation, plot_resp, plot_importance),
  width = 10,
  height = 9,
  dpi = 120
)

ggsave(
  'output/appendix_n2.jpg',
  plot_combination,
  width = 10,
  height = 9,
  dpi = 120
)

################
## Appendix M ##
################

####
## Set up template ---
####

pd <- position_dodge(0.75)

plot_canvas_typedem <- function(df=df,x=name,y=estimate,se=se, group = group){
  ggplot(data = df, aes(y = estimate, x = name, group = group, colour = group)) + 
    geom_point(position = pd) + 
    geom_errorbar(aes(ymin = estimate-1.96*se, ymax = estimate+1.96*se),
                  width = 0,
                  position = pd,
                  size = 1) +
    coord_flip() +
    ylim(-0.2,0.5) +
    labs(x = "", 
         y = "Marginal Effects",
         color = "black") +
    geom_hline(yintercept = 0,
               linetype = "dashed") +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
          legend.position = "bottom",plot.title = element_text(size=14),
          axis.text.y = element_text(colour = "black",size=14),
          axis.text.x = element_text(colour = "black",size=14),
          text=element_text(family="Times",size=14,color="black")) +
    scale_colour_manual(
      name = "", 
      breaks = c("Polity IV","Bjørnskov and Rode"),
      labels = c("Polity IV","Bjørnskov and Rode"),
      values = c("Grey", "Black"))
}

###
# Polity ---
###

df_aut <- read_rds("morning_after_dataset_individuallevel.rds") %>% filter(lag_polity2 < 7 & leader == 0 & !is.na(purged_y1))

###
# Marginal effects ---
###

# Responsiblity ---

responsibility_felm <-  felm(purged_y1 ~ minister_type/coupattempt_presentyear +
                               gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                               I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility_polity <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                          "Minister of Natural Resources","Other type of minister"),
                                 estimate = c(responsibility_felm$coefficients[10],
                                              responsibility_felm$coefficients[11],
                                              responsibility_felm$coefficients[12],
                                              responsibility_felm$coefficients[13],
                                              responsibility_felm$coefficients[14]),
                                 se = c(responsibility_felm$coefficients[10,2],
                                        responsibility_felm$coefficients[11,2],
                                        responsibility_felm$coefficients[12,2],
                                        responsibility_felm$coefficients[13,2],
                                        responsibility_felm$coefficients[14,2]
                                 )) %>% mutate(name = fct_relevel(name,
                                                                  "Minister of Natural Resources",
                                                                  "Other type of minister",
                                                                  "Minister of Finance",
                                                                  "Minister of Foreign Affairs",
                                                                  "Minister of Defense"),
                                               group = "Polity IV"
                                 )


# Affiliation ---

df_party <- df_aut %>% filter(party_group != "unknown") %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm <- felm(purged_y1 ~ party_group/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation_polity <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(affiliation_felm$coefficients[8],
                                        affiliation_felm$coefficients[9],
                                        affiliation_felm$coefficients[10]),
                           se = c(affiliation_felm$coefficients[8,2],
                                  affiliation_felm$coefficients[9,2],
                                  affiliation_felm$coefficients[10,2]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"),
             group = "Polity IV"
)


# Importance ---

df_importance <- df_aut %>% filter(importance_factor != "0")

importance_felm <- felm(purged_y1 ~ importance_factor/coupattempt_presentyear +
                          gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                          I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance_polity <- tibble(name = c("Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                   "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                          estimate = c(importance_felm$coefficients[10],
                                       importance_felm$coefficients[11],
                                       importance_felm$coefficients[12],
                                       importance_felm$coefficients[13],
                                       importance_felm$coefficients[14]),
                          se = c(importance_felm$coefficients[10],
                                 importance_felm$coefficients[11,2],
                                 importance_felm$coefficients[12,2],
                                 importance_felm$coefficients[13,2],
                                 importance_felm$coefficients[14,2]
                          )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
             group = "Polity IV"
)


# Experience ---

df_experience <- df_aut %>% filter(importance_factor != "0.5")

quantiles_experience <- df_aut %>% summarise(x = quantile(experience, c(0.25, 0.5, 0.75)), experience = c(0.25, 0.5, 0.75))

experience_group <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_polity <- tibble(name = c("Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years"),
                                estimate = c(experience_group$coefficients[6],
                                             experience_group$coefficients[7],
                                             experience_group$coefficients[8],
                                             experience_group$coefficients[9]),
                                se = c(experience_group$coefficients[6,2],
                                       experience_group$coefficients[7,2],
                                       experience_group$coefficients[8,2],
                                       experience_group$coefficients[9,2]
                                )) %>% mutate(name = fct_relevel(name,
                                                                 "Over 6 years","4-6 years",
                                                                 "2-3 years",
                                                                 "Less than 2 years"),
                                              group = "Polity IV"
                                )


# Combination ---

medianbyyear <- df_experience %>% group_by(country_isocode,year) %>% summarize(experience_median = median(experience,na.rm=TRUE))

df_combination <- df_experience %>% left_join(.,medianbyyear,by=c("country_isocode","year")) %>% filter(party_group != "unknown") %>%
                                                                                                  mutate(experience_bin = ifelse(experience >= experience_median, 1,0),
                                                                                                        important_bin = if_else(importance_numeric %in% c(3,4),1,0),
                                                                                                        loyal_bin = if_else(experience_bin == 1 & party_group == "leader_party",1,0),
                                                                                                        typology = case_when(important_bin == 0 & loyal_bin == 0 ~ 1,
                                                                                                                             important_bin == 1 & loyal_bin == 0 ~ 2,
                                                                                                                             important_bin == 0 & loyal_bin == 1 ~ 3,
                                                                                                                             important_bin == 1 & loyal_bin == 1 ~ 4
                                                                                                        )
)

combination_felm <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination_polity <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                    "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                           estimate = c(combination_felm$coefficients[6],
                                        combination_felm$coefficients[7],
                                        combination_felm$coefficients[8],
                                        combination_felm$coefficients[9]),
                           se = c(combination_felm$coefficients[6,2],
                                  combination_felm$coefficients[7,2],
                                  combination_felm$coefficients[8,2],
                                  combination_felm$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                            "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"),
                                         group = "Polity IV"
                           )


# Get N's ---

n_exp <-  df_experience %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_presentyear) & !is.na(minister_type)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

####
# Bjørnskov and Rode ---
####

df_aut <- read_rds("morning_after_dataset_individuallevel.rds") %>% filter(democracy_lag_br == "Autocracy" & leader == 0 & !is.na(purged_y1))

###
# Marginal effects ---
###

# Responsiblity ---

responsibility_felm <-  felm(purged_y1 ~ minister_type/coupattempt_presentyear +
                               gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                               I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_aut) %>% summary()

coef_responsibility_dd <- tibble(name = c("Minister of Defense","Minister of Finance","Minister of Foreign Affairs",
                                       "Minister of Natural Resources","Other type of minister"),
                              estimate = c(responsibility_felm$coefficients[10],
                                           responsibility_felm$coefficients[11],
                                           responsibility_felm$coefficients[12],
                                           responsibility_felm$coefficients[13],
                                           responsibility_felm$coefficients[14]),
                              se = c(responsibility_felm$coefficients[10,2],
                                     responsibility_felm$coefficients[11,2],
                                     responsibility_felm$coefficients[12,2],
                                     responsibility_felm$coefficients[13,2],
                                     responsibility_felm$coefficients[14,2]
                              )) %>% mutate(name = fct_relevel(name,
                                                               "Minister of Natural Resources",
                                                               "Other type of minister",
                                                               "Minister of Finance",
                                                               "Minister of Foreign Affairs",
                                                               "Minister of Defense"),
                                            group = "Bjørnskov and Rode"
                                            
                              )


coef_responsibility_both <- rbind(coef_responsibility_polity,coef_responsibility_dd)

plot_resp_both <- plot_canvas_typedem(coef_responsibility_both) + ggtitle("H3a: Ministerial responsibility")

# Affiliation ---

df_party <- df_aut %>% filter(party_group != "unknown") %>% mutate(party_group = as.factor(recode(party_group,"independent"="No party affiliation","leader_party"="From the leader's party","other_party"="From another party")))

affiliation_felm <- felm(purged_y1 ~ party_group/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                           I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_party) %>% summary()

coef_affiliation_dd <- tibble(name = c("From another party","From the leader's party","No party affiliation"),
                           estimate = c(affiliation_felm$coefficients[8],
                                        affiliation_felm$coefficients[9],
                                        affiliation_felm$coefficients[10]),
                           se = c(affiliation_felm$coefficients[8,2],
                                  affiliation_felm$coefficients[9,2],
                                  affiliation_felm$coefficients[10,2]
                           )
) %>% mutate(name = fct_relevel(name,
                                "From the leader's party",
                                "No party affiliation",
                                "From another party"),
             group = "Bjørnskov and Rode"
             
)

coef_affiliation_both <- rbind(coef_affiliation_polity,coef_affiliation_dd)

plot_affiliation_both <- plot_canvas_typedem(coef_affiliation_both) + ggtitle("H2b: Affiliation")

# Importance ---

df_importance <- df_aut %>% filter(importance_factor != "0")

importance_felm <- felm(purged_y1 ~ importance_factor/coupattempt_presentyear +
                          gdp_cap_pwt_ln + pop_pwt_ln + experience + I(experience^2) +
                          I(experience^3) | country_isocode+year | 0 | country_isocode , data=df_importance) %>% summary()

coef_importance_dd <- tibble(name = c("Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                   "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
                          estimate = c(importance_felm$coefficients[10],
                                       importance_felm$coefficients[11],
                                       importance_felm$coefficients[12],
                                       importance_felm$coefficients[13],
                                       importance_felm$coefficients[14]),
                          se = c(importance_felm$coefficients[10],
                                 importance_felm$coefficients[11,2],
                                 importance_felm$coefficients[12,2],
                                 importance_felm$coefficients[13,2],
                                 importance_felm$coefficients[14,2]
                          )
) %>% mutate(name = fct_relevel(name,
                                "Junior minister\nor other low-\nranking post","Low-ranking\nminister","Medium-ranking\nminister",
                                "Vice president,\ndeputy prime minister,\ntop minister","Prime minister/President\n(not leader)"),
             group = "Bjørnskov and Rode"
             
)

coef_importance_both <- rbind(coef_importance_polity,coef_importance_dd)

plot_importance_both <- plot_canvas_typedem(coef_importance_both) + ggtitle("H3b: Importance in cabinet") + ylim(-0.1,0.7)

# Experience ---

df_experience <- df_aut %>% filter(importance_factor != "0.5")

quantiles_experience <- df_aut %>% summarise(x = quantile(experience, c(0.25, 0.5, 0.75)), experience = c(0.25, 0.5, 0.75))

experience_group <- felm(purged_y1 ~ as.factor(experience_group)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode, data=df_experience) %>% summary()

coef_experience_dd <- tibble(name = c("Less than 2 years","2-3 years",
                                         "4-6 years","Over 6 years"),
                                estimate = c(experience_group$coefficients[6],
                                             experience_group$coefficients[7],
                                             experience_group$coefficients[8],
                                             experience_group$coefficients[9]),
                                se = c(experience_group$coefficients[6,2],
                                       experience_group$coefficients[7,2],
                                       experience_group$coefficients[8,2],
                                       experience_group$coefficients[9,2]
                                )) %>% mutate(name = fct_relevel(name,
                                                                 "Over 6 years","4-6 years",
                                                                 "2-3 years",
                                                                 "Less than 2 years"),
                                              group = "Bjørnskov and Rode"
                                              
                                )

coef_experience_both <- rbind(coef_experience_polity,coef_experience_dd)

plot_expgroup_both <- plot_canvas_typedem(coef_experience_both) + ggtitle("H2a: Experience")

# Combination ---

medianbyyear <- df_experience %>% group_by(country_isocode,year) %>% summarize(experience_median = median(experience,na.rm=TRUE))

df_combination <- df_experience %>% left_join(.,medianbyyear,by=c("country_isocode","year")) %>% filter(party_group != "unknown") %>%
                                                                                                 mutate(experience_bin = ifelse(experience >= experience_median, 1,0),
                                                                                                        important_bin = if_else(importance_numeric %in% c(3,4),1,0),
                                                                                                        loyal_bin = if_else(experience_bin == 1 & party_group == "leader_party",1,0),
                                                                                                        typology = case_when(important_bin == 0 & loyal_bin == 0 ~ 1,
                                                                                                                             important_bin == 1 & loyal_bin == 0 ~ 2,
                                                                                                                             important_bin == 0 & loyal_bin == 1 ~ 3,
                                                                                                                             important_bin == 1 & loyal_bin == 1 ~ 4
                                                                                                        )
)

combination_felm <- felm(purged_y1 ~ as.factor(typology)/coupattempt_presentyear +
                           gdp_cap_pwt_ln + pop_pwt_ln | country_isocode+year | 0 | country_isocode , data=df_combination) %>% summary()

coef_combination_dd <- tibble(name = c("Low responsibility and weak signs of loyalty","High responsibility and weak signs of loyalty",
                                    "Low responsibility and strong signs of loyalty","High responsibility and strong signs of loyalty"),
                           estimate = c(combination_felm$coefficients[6],
                                        combination_felm$coefficients[7],
                                        combination_felm$coefficients[8],
                                        combination_felm$coefficients[9]),
                           se = c(combination_felm$coefficients[6,2],
                                  combination_felm$coefficients[7,2],
                                  combination_felm$coefficients[8,2],
                                  combination_felm$coefficients[9,2]
                           )) %>% mutate(name = fct_relevel(name,
                                                            "High responsibility and strong signs of loyalty","Low responsibility and strong signs of loyalty",
                                                            "High responsibility and weak signs of loyalty","Low responsibility and weak signs of loyalty"),
                                         group = "Bjørnskov and Rode"
                                                            
                           )



coef_combination_both <- rbind(coef_combination_polity,coef_combination_dd)

plot_combination_both <- plot_canvas_typedem(coef_combination_both) + ggtitle("H4: Combinations of traits")

# Get N's ---

n_exp <-  df_experience %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_affliation <-  df_party %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_resp <-  df_aut %>% filter(!is.na(coupattempt_presentyear) & !is.na(minister_type)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_importance <-  df_importance %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))

n_comb <-  df_combination %>% filter(!is.na(coupattempt_presentyear)) %>% group_by(coupattempt_presentyear) %>% summarize(n = as.numeric(n())) %>% ungroup() %>%
  mutate(sum = sum(n))


###
# Print ---
###

ggsave(
  'output/appendix_m1.jpg',
  gridExtra::grid.arrange(plot_expgroup_both, plot_affiliation_both, plot_resp_both, plot_importance_both),
  width = 11,
  height = 10,
  dpi = 120
)

ggsave(
  'output/appendix_m2.jpg',
  plot_combination_both,
  width = 11,
  height = 10,
  dpi = 120
)